Vortex signatures in annular Bose-Einstein condensates 
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We consider a Bose-Einstein condensate confined in a "Mexican hat" potential, with a quartic 
minus quadratic radial dependence. We find conditions under which the ground state is annular 
in shape, with a hole in the center of the condensate. Rotation leads to the appearance of stable 
multiply-quantized vortices, giving rise to a superfiuid fiow around the ring. The collective modes 
of the system are explored both numerically and analytically using the Gross-Pitaevskii and hydro- 
dynamic equations. Potential experimental schemes to detect vorticity are proposed and evaluated, 
which include measuring the splitting of collective mode frequencies, observing expansion following 
release from the trap, and probing the momentum distribution of the condensate. 

PACS numbers: 03.75.Lm, 03.75.Kk, 67.40.Vs 



I. INTRODUCTION 

The 1995 discovery of Bose-Einstein condensation in 
dilute ultracold gases P, IH prompted a flurry of activ- 
ity investigating these fascinating quantum systems. In 
the experiments, ultracold atoms are confined by mag- 
netic or optical forces, which typically can be represented 
theoretically by a harmonic potential. Lately, however, 
an increasing amount of interest has focused on Bose con- 
densates in more exotic trapping potentials. For exam- 
ple, Gupta et al. Q recently produced a condensate in 
a ring-shaped magnetic waveguide, while experiments at 
Ecole Normale Superieurc 5, 6] have been conducted in 
overlapping magnetic and optical dipole traps, leading 
to a potential which can be approximated by a function 
which is quadratic plus quartic in the radial coordinate. 
This trap potential is especially of interest when the con- 
densate is rotating, since, in contrast to the harmonic 
case, one can in principle achieve arbitrarily high angular 
velocities. This leads to interesting new configurations, 
for example a vortex lattice with a hole at its center, and 
for even higher angular velocities, annular condensates 
containing a multiply-quantized "macrovortex" . These 
equilibrium configurations have been studied in a num- 
ber of papers 0,13,13' where the angular velocities for 
transitions between the different states, defining a phase 
diagram, were calculated in Refs. 0, 0, lia^lTsL IT3 |. 
In addition to the static properties, the dynamics have 
also been addressed, with the collective mode excitations 
considered in Refs. [l3.[l5j|. 

Alternatively, one could consider the case where, in ad- 
dition to the quartic term, there is a negative harmonic 
term. Then the radial trap potential has the form of a 
"Mexican hat" , and the condensate can be annular even 
when there is no rotation. Similar trap configurations 
have been studied elsewhere [l^ d, ll^l ■ Of particu- 
lar interest, however, is the behavior of the system under 
rotation, since vortices can enter the condensate and form 
a stable, multiply-quantized vortex at the center. For 
example, an important question raised in this configura- 



tion concerns the stability of the vortex in the presence 
of a non-rotating thermal cloud at finite temperatures — 
does the vortex decay by spiraling towards the edge as 
in a harmonic trap |23, |21| , or is it metastable? This is 
related to the well-known problem of stability of persis- 
tent currents in superfluids [2^ , and so is of interest 
from a fundamental perspective. 

In order to tackle these types of questions experimen- 
tally it is important to possess some way of detecting 
the presence of a vortex. In a harmonically-trapped con- 
densate this is relatively straightforward, since the fluid 
circulation around the vortex core creates a hole in the 
density, which can be detected by standard absorption 
imaging. Since the optical resolution is usually not suffi- 
cient to detect vortices in situ, a period of free expansion 
following the release of the condensate from the trap is 
generally required. However, in a trapped annular con- 
densate vortices reside in the center where the density is 
very low irrespective of the circulation. Hence, imaging 
of a vortex in situ is precluded, and it is not even obvious 
a priori whether it is possible to resolve the vortex after 
expansion. A more detailed analysis is thus required in 
order to evaluate the feasibility of this and other possible 
schemes for vortex detection. 

In this paper we shall study two-dimensional con- 
densates, where the trap potential has a quartic minus 
quadratic radial dependence. Most of our results for thin 
annular condensates, however, also apply to more general 
Mexican hat potentials. Two dimensions in this context 
corresponds to either a cylindrical condensate, where the 
axial length scale greatly exceeds the radial size, or a 
thin disk-shaped condensate where the system is tightly 
trapped in the axial direction such that dynamics along 
this coordinate are "frozen out" and can be neglected. 
Both geometries are of experimental interest. We also 
focus exclusively on condensates in the mean field regime 
at zero temperature, where the properties can be accu- 
rately described by the Gross-Pitaevskii (GP) equation. 
Nevertheless, we should mention that the possibility of 
obtaining a thin annulus, essentially a one-dimensional 
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system, opens up other interesting regimes beyond mean 
field, such as a quasi-condensate or Tonks- Girardeau gas 

We begin by studying the equilibrium states, where 
both numerical solutions and analytical results are pre- 
sented. We then characterize the lowest energy excita- 
tions, or collective modes, of the condensate. In partic- 
ular, we use a hydrodynamic model to derive analyti- 
cal expressions in the limit of thin annuli, as well as to 
find the mode frequencies numerically. We also compare 
these results to numerical simulations of the GP equa- 
tion. This allows us to discuss the use of collective modes 
to detect vorticity; namely the splitting between counter- 
propagating surface and acoustic modes of the conden- 
sate. We then consider the expansion of the condensate 
after switching off the confining potential, finding that 
the hole due to the vortex is recovered in the expanded 
density, providing a simple means to detect the presence 
of the vortex in experiments. Finally, we discuss the mo- 
mentum distribution of an annular condensate, finding 
that in the limit of thin annuli the function approaches 
a simple analytical form. 



II. EQUILIBRIUM CONFIGURATIONS 

For a dilute Bose-Einstein condensate in the limit 
of zero temperature, the time-dependent evolution of 
the condensate wavefunction \E'(r,t) is described by the 
Gross-Pitaevskii (GP) equation 
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The interactions between atoms are represented by the 
parameter (73D = AnNfi^a/M , for N atoms of mass M, 
with s-wave scattering length a. Note that, in addition 
to the dynamics of the condensate, the time-independent 
stationary states can be found by substituting \E'(r,t) — 
'^[r)e~'^^*^^ into Eq. whereupon the left-hand side 
becomes /j.'ijj with /x the chemical potential. The external 
potential in Eq. ^ is of the form 



where the evolution in the axial coordinate is neglected. 
Physically this can correspond to the case of a cylindri- 
cal condensate, where ujz ^ The interaction param- 
eter then becomes, g = AirN'a, in harmonic oscillator 
units, where N' = N/Z represents the number of atoms 
per unit length along the axial direction. An alternative 
two-dimensional configuration involves tight axial con- 
finement, such that hujz 3> 53_d|^'(0)P in physical units. 
In this geometry, which has been realized experimentally 
[25I I , the axial dynamics are "frozen out" . The axial de- 
pendence of the wavefunction is then gaussian, and the 
2D GP equation ^ can be used with N' = N/{V2Trdz) 

M 

We first analyze a condensate without a vortex using 
the Thomas-Fermi (TF) approximation, where the spa- 
tial derivative of 5* in Eq. Q is neglected, so that the 
density is given by n{r) = |5'(r)p — {fi — Vc^t)/g for 
n{r) > 0, and n{r) — otherwise. The radii of the con- 
densate are evaluated by finding where the density goes 
to zero. If /i > then there is a single solution of the 
resulting quadratic equation, i?^ = (1-1- VI + 8A/x)/(2A), 
which corresponds to a condensate without a hole. On 
the other hand, if /i < there are two solutions i?i 2 = 
(l=F-\/l + 8A/x)/(2A), and the condensate is annular with 
a hole in the center. In this case one can write the TF 
condensate density as 



n{r) = ^{Rl-r'){r'-Rl)., 
2.9 



(4) 



for n{r) > 0. Calculating the radii and chemical poten- 
tial more explicitly requires the condition that the wave- 
function is normalized to unity Jdrn{r) — 1. The re- 
sulting radii can be conveniently expressed as the sum 
and difference of the squares = i?| ± i?^ , yielding 
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where we have defined a new parameter ?/ in Eq. ((HJ . The 
chemical potential is then 
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where di = y/H/{MuJi), with lu±_ and ujz the radial and 
axial trap frequencies respectively. In the following we 
shall mostly use harmonic oscillator units where distance. 



time and energy are expressed in units of dj 



and 



hujj_ respectively. We also take a — —1, where the case of 
(7 = 1 has already been extensively studied in a number 
of papers 0, i, Ulli IH El El El El . 

We explicitly consider the GP equation in 2D 
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This yields the condition 77 < 1 for an annular conden- 
sate, which we will always assume in the following. The 
limit 77 — > corresponds to the case where the width of 
the annulus becomes much smaller than the radius. 

The range of validity of the TF approximation for an 
annular condensate is defined by the condition ^ d 
figj. where ^ — {2gnjnax)~^^^ is the healing length at the 
density maximum, while d = R2 ~ Ri is the width of the 
annulus. Using the above expressions yields the explicit 
relations ^ = 2^/\/r] and d = + rj — ^\ — r\) j V2A. 

The condensate can also contain a single centered vor- 
tex with circulation F = 2TTvh/M , where the well-known 
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requirement of quantized circulation in superfluids |27| 
constrains the vorticity to integer values ly = 0,1, 2, ... . 
The equilibrium condensate density is represented within 
the TF approximation as 



where we have defined the second-order term separately 
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gn{r) = ^ - i I ^ - + Ar 



(8) 



where fi is the chemical potential in the laboratory (non- 
rotating) frame and the jr^ term corresponds to the 
centrifugal potential barrier resulting from the fluid flow 
around the vortex. Performing an analysis similar to 
that of the non-rotating condensate gives the following 
relations 
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Interestingly, this term can also be derived by calculat- 
ing the kinetic energy from the circulation of the vortex 
J dr nv'^ / (2r'^) , where n is approximated by the density 
without a vortex Q). 

The energy in a frame rotating with angular veloc- 
ity f2 is related to that in the laboratory frame by 
E'{i') = £{1/) — fii/. For a particular there exists a 
vorticity v that minimizes this energy, identifying the 
ground state for each rotation rate. Treating j/ as a con- 
tinuous quantity (valid for i> 1) and using Eq. H17|) to 
second order gives the following condition for the ground 
state vorticity 
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Approximate analytical solutions of Eqs. H9I11|I can be 
obtained using the expansion XR\ = 1-1- X^n ^"(Ai^)^", 
which is most appropriate for Ai^ ^ 1. As we shall see 
this regime is the most interesting for our purposes, since 
it corresponds to the region where a multiply-quantized 
vortex is stable. Substituting into Eqs. H9I11I) yields the 
following results 
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Another important quantity is the energy of the con- 
densate in the presence of a vortex. This is represented 
in terms of an energy per particle as an integral 



E / dr 



-V^ - + Xr^ + gn) ^ . (15) 



Within the TF approximation the energy becomes 

72 



E =^ dr - r'^ + Ar^* + gn ] n . (16) 



Eq. H16|l can be approximated analytically using the TF 
density ||SJ) and the expansions I|12I14|I . This eventually 
yields 
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+ ^rj^ + E,{Xvf + 0[{Xvf], (17) 



In the limit rj —t Q this expression becomes particularly 
simple, since Ei ~ 1, so that Xv ~ Vl/2. This can also be 
easily derived by equating the velocity due to the rotation 
at the radius of the annulus, v — QR, to the fluid velocity 
due to the vortex, v = v/R, where R^ = R\/2 = 1/(2A). 

These analytical results can be compared to solutions 
of the GP equation in the rotating frame 



1 . 



Xr^ 



(20) 

where = i{yd/dx — xd/dy) is the axial component of 
the angular momentum operator. To find the equilibrium 
configuration for each Q, we solve Eq. H20|l numerically by 
propagating in imaginary time (i.e. making the replace- 
ment t — > ^it), starting from the il = ground state. 
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FIG. 1: Density profiles of a condensate, found by numerical 
solution of the GP equation with g = 1000 and A = 0.01 
{g = 0.726) for (a) n = 0.2 (in units of tj^), (b) Q. = 0.4, and 
(c) Q. — 0.5; (d)-(f) show the corresponding phase profiles for 
each Q. The length scale of each box is 28 x 28 (in units of 
d±). 
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This technique yields the configurations shown in 
Fig. ^ One can write the condensate wavefunction as 
\I'(r) = ^/rl(r) exp[iS{r)], where n{r) = |\['(r)|^ is the 
density and S{r) is the phase. Fig. ^ shows the density 
and phase at three different rotation rates for a particular 
set of parameters {g = 1000, A = 0.01, t] = 0.726). For 

— 0.2 (Fig.^a)) the equilibrium state consists of a sta- 
ble multiply-quantized vortex with = 11, as can be seen 
from plotting the phase (Fig. ^d)) where, by traversing 
the center in a closed loop, one finds 11 cycles between 
5 = (black) and S ^ 2tt (white). Atn = 0.4 (Figs.mb) 
and (e)) the multiply-quantized vortex is again stable, 
but with v = 22. However, at = 0.5 (Figs, ^^c) and 
(f)) a multiply- quantized vortex at the center of the con- 
densate is surrounded by two rings of singly-quantized 
vortices. The outer ring is visible in the density plot since 
the vortices are embedded in the condensate. This tran- 
sition between a "macrovortex" and a "lattice with hole" 
state was previously found in Ref. JJJ for 3D configu- 
rations. This is similar to the transition in quartic plus 
quadratic traps discussed in Refs. [ifl llll 13 ; ex- 

cept that in the latter case the transition takes place at 
a much higher frequency and in the opposite direction. 
To find the critical rotation rate for the transition in the 
present configuration would require a detailed analysis, 
possibly using techniques similar to Refs. [T3|[l3|. 

In the remainder of the paper we shall focus on the 
macrovortex state, (i.e. one where there is a single, 
multiply-quantized vortex in the center of the hole) where 
in all cases we have checked that the macrovortex is sta- 
ble in the rotating frame by numerically solving the GP 
equation. The analytical results derived above can then 
be used, and compared to numerical solutions of the GP 
equation. In order to facilitate such a comparison, it 
is simpler to solve the GP equation in the non-rotating 
frame with an effective potential Vcy:t{r, z = 0)+i'^ / (2r^). 
The energy can then be calculated from the density us- 
ing Eq. (|15|l and is compared to the analytical expres- 
sion (|17|) in Fig. [3 One sees good agreement for v < 20 
(corresponding to Xv < 0.2), which is the region where 
Eq. ((T7|) is expected to be vahd. To find the ground state 
for a particular rotation rate f2 one can, as before, simply 
identify the v that minimizes the rotating frame energy 
E'{i') = E{v) — vVt. The resulting v vs. plot is com- 
pared to the analytical result H19|) in the inset of Fig. |21 
again showing good agreement for small Xv. 

Before continuing, it is interesting to relate our work 
to more general "Mexican hat" type potentials — for ex- 
ample, a harmonic plus gaussian potential [T^ which ex- 
perimentally would correspond to magnetic trap overlaid 
with a tightly-focussed blue-detuned laser beam. In par- 
ticular, for a thin annulus, when the condensate is con- 
fined close to the potential minimum at radius i?, the 
radial trap potential can be expanded to second order in 
r — i? as 

Vir)^^u;Ur~Rr + 0[ir-Rf] (21) 
where we have neglected a constant term which does 
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FIG. 2: Energy (in units of huj±) of the condensate with a 
central vortex as a function of vorticity v for the same pa- 
rameters as Fig.0 The circles plot numerical solutions of the 
GP equation, while the solid line is the analytical estimate 
l)17|l . Inset: vorticity of the ground state in a frame rotating 
with angular velocity Q (in units of uj±). Circles are numerical 
GP solutions, while the solid line plots the analytical result 



not affect the dynamics. For our potential Q we find 
LUr ~ \/2ujj_ and R = dj_/V2X (in physical units). In 
the following, the results presented for the thin annulus 
limit 77 ^ 1 are therefore applicable to more general po- 
tentials, where finding ujr and R from the corresponding 
expansion H21|) allows a precise mapping onto our system. 



III. COLLECTIVE MODES 

In the previous section it was found that, within the TF 
approximation, a hole is present within a non-rotating 
condensate for 77 < 1. Considering a condensate rotat- 
ing with sufficiently small angular velocities Jl, we found 
equilibrium configurations where multiply-quantized vor- 
tices are stable. An interesting question concerns how 
these vortices could be detected experimentally, since the 
density profiles of a condensate with no vortex or small 
values of i> are very similar, making it very difficult to dis- 
tinguish between the different states on this basis alone. 
One possibility is to study the collective modes of the 
system, since differences in the vorticity may be man- 
ifested in the frequencies, in a similar way to those in 
a harmonic trap [23, 123, yOl ■ The problem of collective 
modes in annular condensates is also interesting in its 
own right, as similar investigations in a rapidly rotating 
quadratic plus quartic trap have shown [l5| . This section 
will therefore detail similar analyses for a quartic minus 
quadratic confinement, with the specific application to 
vortex detection discussed at the end. 

One approach to finding collective mode frequencies is 
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to rewrite the GP equation in the rotating frame H20() 
in terms of the density and phase, and use the TF ap- 
proximation to obtain hydrodynamic equations for these 
quantities. Then hnearizing for small fluctuations about 
equilibrium, SS, 5n, yields 



d_ 
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6n - 
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(22) 
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where uq is the TF equilibriimr density from Eq. 
For density and phase modulations of the form Sn^ SS cx 
gi(m0-wt) ^Yiis gives 6n = i(tj + mQ — miy/r^)SS and 
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rgriQ—SS 



0, (24) 



where co is the frequency of the mode in the rotating 
frame. 

It is instructive to first consider the case of a non- 
rotating condensate where = and v = 0. Simi- 
larly to Ref. 15], one can introduce the new variable 
( = [r^ - R\/2)/{R^_/2), so that the equilibrium den- 
sity Q becomes gn^ = Ai?i(l - C^)/8 and -1 < C < 1- 
Substituting into Eq. H24|) then gives 



1 + 
d 



SS 



d 



{l + rjOil-Og^SS 



0. (25) 



Note that this equation, and therefore the mode fre- 
quency, depends only upon the parameter rj and the 
square of the azimuthal quantum number to, the posi- 
tive and negative to modes being degenerate for symme- 
try reasons. In the limit of rj —^ (i.e. when the width 
of the annulus is small compared to its radius, as evident 
from the relation t] = R^/R^) the equation becomes 



LO^SS^ 



d_ 



(i-C^)|^5 



0, 



(26) 



which is Legendre's equation with eigenfunctions given 
by Legendre polynomials Pj{C) s-^d eigenvalues = 
j{j + 1), where j = 0, 1, 2, . . . is the number of radial 
nodes ^V\. Hence, different to states are degenerate in 
this limit 32]. The same result for the frequencies can 
also be derived starting from the expanded potential H21() 
with the corresponding one-dimensional hydrodynamic 
equation for the radial coordinate. This implies that the 
modes can be pictured as consisting of different parts of 
the annular condensate undergoing radial oscillations in 
a one-dimensional harmonic well, providing a simple ex- 
planation for the mode frequencies and the degeneracy 
of different to states. 



Apart from these analytical considerations, one can 
also solve Eq. numerically using a shooting method 
|33| to find the mode frequencies ui and functions SS{(). 
Fig. inja) shows the high- lying (j = 1) m = 0, 1, 2 mode 
frequencies as a function of rj for the non-rotating {Q = 0, 

= 0) condensate in the presence of the hole {i] < 1). 
One sees that in the limit of ?7 ^ all of the mode fre- 
quencies tend to w = \/2, consistent with the above anal- 
ysis. With increasing 77, however, the frequencies deviate 
from this value, with a particularly marked decrease in 
the TO = mode as one approaches rj — 1. In addition, it 
becomes increasingly difficult to find numerical solutions 
near to 77 = 1. Numerical results from simulating the 
time-dependent GP equation © for each mode are also 
presented in Fig.|2Ia). The GP results closely follow the 
hydrodynamic equations up to 77 = 0.7, above which the 
GP results start to deviate significantly. This discrep- 
ancy is related to the fact that when the inner radius 
approaches zero at 77 — > 1 the contrast between the den- 
sity profiles for the TF approximation (where the density 
goes to zero suddenly) and from solving the GP equation 
(where the density tails off more gradually) becomes in- 
creasingly significant, leading to differences in the nature 
of the modes near to the center. In addition, the transi- 
tion between hole and non-hole states in the GP solution 
is more gradual, avoiding the singular behavior in the 
hydrodynamic equation at rj = 1. 
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FIG. 3: (a) Frequency (in units of uj±) of the high-lying 
(j = 1) 771 = (solid line) m = 1 (dashed) and m = 2 
(dotted) modes as a function of rj for a non-rotating conden- 
sate [Q. — 0), from numerical solution of the hydrodynamic 
equation II25II . Also plotted are numerical solutions of the 
Gross-Pitaevskii equation ^ ai g — 1000 for m — {o), 
m — 1 (□) and m = 2 (A), (b) Frequencies of the low-lying 
[j = Q) m = 1 and m = 2 modes, where the grey lines show 
the analytical result 1281 for each m, while the other lines and 
points are labelled similarly to (a). 



6 



The low-lying (j = 0) m > modes can also be inves- 
tigated, and the m — 1 and m = 2 modes frequencies are 
plotted as a function of rj in Fig. I^Ib). Comparison be- 
tween the hydrodynamic and GP results again show very 
good agreement. Similar modes were found in a rapidly 
rotating quadratic plus quartic traps in Ref. 15], and in 
the thin annulus limit rj ^ 1 correspond to compressional 
modes directed azimuthally around the annulus, in con- 
trast to the high-lying modes which correspond to shape 
oscillations of the annulus. This physical explanation is 
the basis of a simple analysis, starting from the relation 
Lu = cq, where for a narrow annulus q — m/R, with the 
mean radius R ~ i?+ / ^/2. The speed of sound, c, is given 
by the relation — nidfj,/dni, where ni = N'/{2tt) is 
the density of the condensate integrated over the radial 
direction. Rewriting in terms of ry yields the equation 



0J^ = 



2m 



3\Rl''dri 



(27) 



Substituting Eqs. (0) and for a non-rotating conden- 
sate, Ai/ = 0, gives c — ri/^/l2\ and 



1^ 



(28) 



The mode frequencies for m = 1 and to = 2 found from 
Eq. (|28|l are plotted as a function of 77 in Fig. Elb). This 
analytical estimate is exact for thin annuli ,32.] . and in- 
deed the comparison to the hydrodynamical results shows 
excellent agreement for ry < 0.4. 

So far we have only studied the non-rotating case where 
il = 0. In a frame rotating with a non-zero angular fre- 
quency, 57 ^ 0, we demonstrated earlier that for suffi- 
ciently small f2 the ground state consists of a multiply 
quantized vortex v > 1. These stationary states can be 
described analytically within the TF approximation us- 
ing Eq. ||SJ| for the density. From the normalization con- 
dition for the density one can find the chemical potential 
and energy for each v, followed by the v state which min- 
imizes the energy in a frame rotating with frequency Q. 
Substituting the resulting v and no into Eq. (|24ll and 
solving numerically then yields mode frequencies for this 
particular ft. 

We find that the rotating-frame frequencies for the 
high- lying {j — 1) modes change significantly with in- 
creasing n, with a rather large splitting between to — 
+ \m\ and m = —\m\ modes for a given \m\ > 0. In con- 
trast, the dependence of the low- lying (j = 0) mode 
frequencies is less pronounced, and is almost independent 
of the sign of m. In the laboratory frame, however, one 
has a frequency splitting between low-lying modes of op- 
posite TO. Using a sum rule analysis similar to that of 
Ref. 28] gives for this splitting 



Alj = 4\m\\iy[l + 0{T]^)] 



(29) 



where a similar expression was derived using perturba- 
tion theory in Ref. (18j for a thin toroid in a quadratic 



plus gaussian potential. For the thin annulus (77 ^ 1) one 
can use Eqs. (^5)) and to derive a simple equation for 
the frequencies of the low-lying modes in the laboratory 

frame 



UJ± 



V6 



±2Xiy 



(30) 



Note that these expressions can also be derived by us- 
ing the quasi-degeneracy between modes in the rotating 
frame, together with the transformation tJiab = ujj-ot+mQ 
and Eq. (O with Ei~l. 

In Fig. ^a) we plot the laboratory-frame frequencies, 
I'^iabl, of the TO = ±2 low-lying modes as a function of v 
and H, for a thin annulus, where A — 0.00322, g — 680, 
and 77 = 0.3. To place these values into physical con- 
text, note that they would correspond to e.g. 10^ *^Rb 
atoms in a cylindrical condensate 100 /im in length, with 
trap frequency uj± j (27r) ~ 20IIz giving a length unit of 
dj_ — 2A\MD.. This would give an annulus with inner and 
outer radii of 25 \xva. and 34 /im respectively, and is of in- 
terest as it corresponds to feasible parameters in future 
experiments. Eq. (|30|l is also plotted in Fig.^Ja), where 
we see very good agreement with the hydrodynamic re- 
sults for < 40. 
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FIG. 4: (a) Frequencies in the laboratory frame of the low- 
lying m = ±2 modes as a function of v and Vl, for g = 680 
and A = 0.00322. The solid lines corresponds to numerical so- 
lutions of Eq. 1241 , while the dashed lines show the analytical 
expression (1301 for both modes, (b) The quadrupole moment 
{pp" — y^) (normalized to the mean squared radius {x^ -f y^)) 
as a function of time (in units of given by a numerical 

simulation of the GP equation in the laboratory frame. The 
original condensate contains a single v = 1 vortex, and the 
applied perturbation mainly excites the low-lying m = ±2 
modes. 

Above a particular value of v we find that uj of the 
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m = —2 mode becomes negative. From Eq. (I30|l this 
critical value is found to be independent of |m| and given 
by t'crit = ill good agreement with the nu- 

merical hydrodynamic results. This corresponds to the 
point at which the m = — |m| modes become energeti- 
cally unstable, whereupon a source of dissipation (such 
as the presence of a thermal cloud at non-zero temper- 
atures) results in the mode amplitudes increasing with 
time. This can also be seen from the famous Landau 
criterion 26:|, which states that excitations in a super- 
fluid become energetically unstable at flow velocities, f , 
exceeding the speed of sound, c. For a thin annulus sur- 
rounding a vortex the fluid velocity is f = i^/ R, while 
c — rj/y/vlx, immediately leading to the equation for 
Vcrit derived above. 

Another interesting consequence of the fact that the 
mode frequencies pass through zero is related to a pos- 
sible experiment where one applies a small, static per- 
turbation of multipolar form to a condensate containing 
a ~ Vent vortex. This is similar to the experiment 
performed for a vortex lattice in a harmonically-trapped 
condensate [s^. whereupon the negative m modes are 
resonantly excited and the resulting non-linear behavior 
studied. 

An important question relates to the feasibility of ex- 
ploiting the splitting between the positive and nega- 
tive mode frequencies as a diagnostic for the presence 
of a vortex, in analogy with the similar use of the sur- 
face mode splitting in harmonically-trapped condensates 
[2E HOl ■ We test this idea in more detail by perform- 
ing numerical simulations based on the time-dependent 
GP equation. We start with a condensate with the above 
parameters, and with av — 1 vortex in the center. A per- 
turbation is added to the lab. frame trap potential of the 
form 5V = a(x'^ — y^) up to t = 1, after which the po- 
tential is returned to its original form. The subsequent 
oscillations are tracked by calculating various moments 
of the system {x){t) — J dr |>I'(r, t)px(r) as a function 
of time. 

The (a;^ — y^) quadrupole moment (normalized to the 
mean squared radius (a;^ +2/^)) is shown in Fig. 2{b) 
for a = 0.0025. The perturbation predominately excites 
the low-lying m = +2 and m = —2 modes |35|, and 
analysis of the time series shows that modes with (lab. 
frame) frequencies lu = 0.230 and u; — 0.257 are present, 
which agree well with the predictions from the analytical 
estimate (|30|l of uj- — 0.232 and lu+ = 0.258. Impor- 
tantly the splittings also closely agree, within 5 %, allow- 
ing Eq. I|30|) to be used to deduce v from measurements of 
the splitting. Physically the frequency splitting is respon- 
sible for the "beating" between the two modes evident in 
the plot in Fig. ^fb). In experiments these modes would 
be apparent as a density pattern around the annulus, 
consisting of four antinodes with alternate maxima and 
minima. The presence of the splitting would then lead to 
rotation of the position of the antinodes, similar to the 
"precession" of the quadrupole surface modes seen ex- 
perimentally in harmonic traps |29l l30| . The precession 



frequency Wpr is equal to Atj/(2|m|) [26j, corresponding 
to the angular velocity of the reference frame where the 
modes are degenerate. Note that in the annular conden- 
sate the precession frequency is simply 

Wpr = 2Xiyc^n, (31) 

independent of |m|, in contrast to the harmonically 
trapped case where the precession frequency depends 
upon the multipolarity. The second equality is valid for 
large v and small 77, where one can use Eq. I|19|) with 
El — 1. The simulations were repeated for the |m| = 1 
modes, as well as the |m| — 2 modes in the presence of 
a v — 2 vortex, yielding frequencies that agree well with 
the predictions of Eq. l|Sn|l . 

Thus, the splitting of the modes may be observable 
and should match closely with that found from simple 
analytical estimates. However, a limiting factor could be 
the ability of experiments to resolve the density fluctua- 
tions associated with the mode, since this would proba- 
bly be more difficult than observing shape changes as in 
harmonic traps. One could in principle increase the am- 
plitude of the mode by using a large perturbation (large 
a), but as seen in Fig.^J^b) this can lead to nonlinear cou- 
pling of the mode to higher excitations, leading to a rapid 
damping (note that this damping is absent for smaller 
perturbations, e.g. a — 0.0005). So there is a trade off 
involved in choosing the size of the perturbation, and this 
may in practice limit the utility of collective modes for 
vortex detection. 



IV. EXPANSION 

One alternative to using mode frequencies to detect 
vorticity in annular condensates was proposed by Tem- 
pere et al. [T^ . and involves allowing an expanding con- 
densate to overlap with an another expanding condensate 
without a vortex. The resulting interference pattern re- 
veals the presence of a vortex in the first condensate. 
However, this method may be difficult to realize experi- 
mentally, and a much simpler procedure would be to use 
only one condensate and allow it to expand freely after 
switching off the trapping potential. The central ques- 
tion to be addressed in this section is whether imaging 
the resulting density profile would provide a means to 
distinguish between states containing a vortex and those 
without. 

We study this question by starting with the GP equi- 
librium state either with or without a vortex, then propa- 
gating the 2D time-dependent GP equation with t^xt set 
to zero. We study the narrow annulus case mentioned at 
the end of the previous section, where A — 0.00322 and 
•q = 0.3. Fig. [SI shows radial cross-sections through the 
density for = 0, = 1 and u = 2 for various times 
during the expansion. Initially at t = (Fig. (Sfa)) there 
is very little difference between the three configurations, 
illustrating the difficulty in distinguishing between the 
three states in-situ. After release from the trap {t > 0) 




FIG. 5: Density as a function of radial distance through an 
expanding condensate (A = 0.00322, rj — 0.3) with no vortex 
(solid line), a, v — 1 singly-quantized vortex (dashed) and 
a, v — 2 doubly-quantized vortex (dotted) at (a) t — 0, (b) 
t = 4, (c) t = 8, (d) t = 12. The inset in (b) shows the density 
profiles near to the center, illustrating the differences between 
the states in this region. 



we see that the condensate expands radially both out- 
wards and inwards, thereby filhng the annulus. So, even 
at t = 4 (Fig. [SJb)) the different states are, in principle, 
distinguishable due to the density hole produced by the 
centrifugal barrier near to the vortex core. One can also 
see the development of interference fringes away from the 
center. 

Figs. EJc) and (d) show the condensate after further 
expansion. The size of the density dip at the vortex core 
is proportional to the healing length, ^ = 1/ ■^2gn, so that 
the vortex expands along with the condensate due to the 
corresponding decrease in density. The vortex size can be 
quantified in the simulations by measuring the width at 
half the maximum density of the condensate, which for 
1/ = 1 gives radii of 0.75 at t — 8 and 1.2 at t = 12. Using 
the example parameters introduced in the last section, 
this would correspond to vortex core diameters of 3.6 /im 
and 5.6 /xm at 64 ms and 95 ms respectively, approaching 
vortex sizes resolvable in experiments [s^ . 

Note that at later times the i/ = condensate density 
near to the center is relatively uniform, which could pro- 
vide a means to differentiate between different states 
experimentally. For a uniform condensate the radius of 
the vortex core at half-maximum density can be found 
numerically [2^ yielding values of 1.5^ and 3.0^ for i> ~ 1 
and v — 2. This factor of two difference is reproduced 
in our simulations, which at i = 14 yield vortex radii at 
half-maximum of 1.4 and 2.7 for v = I and v ~ 2 respec- 
tively. Using the predicted value ^ = 0.9, from the mea- 
sured peak density, yields good estimates for the vortex 
radii. So in experiments, by measuring the vortex core 
sizes and using simple theoretical relations, one could 




FIG. 6: Density of the expanding condensate (A — 0.00322, 
rj = 0.3) at t = 14 with (a) no vortex (v = 0), (b) a singly- 
quantized vortex {v — 1), and (c) a doubly-quantized vortex 
{u = 2). The length scale of each box is 28 x 28 d±. 



determine whether either a singly or a doubly quantized 
vortex is present |3^ . In general, the larger the vorticity, 
the wider will be the central hole at a given instant after 
release from the trap. However, the quantitative mea- 
surement of large vorticity values could in practice be 
limited by the longer expansion times needed to yield a 
wide enough central plateau, the corresponding reduction 
in density making the experimental imaging difficult. 

Fig.Elshows the 2D density profile at i = 14 for i/ = 0, 
V — 1 and V = 2^ which would be similar to the absorp- 
tion images typically found in experiments. As well as 
showing the density profiles are circularly-symmetric, it 
also illustrates the large differences between the images, 
which adds qualitative support to the idea that it should 
be relatively simple to distinguish between different vor- 
tex states after expansion. 

It is important to note that the central plateau is a 
consequence of repulsive interatomic forces, and for free 
expansion with no interactions a more sharply peaked 
distribution is found. In this latter case the condensate 
wavefunction is given by the analytical expression j26j 



*(r,<) 



\2iiht 



(32) 

where D is the dimensionality, so in this case D = 2. 
Note also that we have returned to using physical units. 
At large times, t ^ MR\/h, the density of the expanded 
condensate is circularly symmetric and tends towards 



(33) 



where |$(fc)p is the momentum distribution of the con- 
densate. This will be discussed in more detail in the next 
section. 



V. MOMENTUM DISTRIBUTION 

A final method for vortex detection involves consider- 
ing the momentum distribution of the condensate, which 
could be probed experimentally by switching off interac- 
tions during expansion (e.g. by exploiting a Feshbach res- 
onance) as discussed in the previous section, or by Bragg 



9 



spectroscopy via the dynamic structure factor 
Mathematically, the momentum distribution |<i>(fc)p is 
obtained by Fourier transforming the equilibrium wave- 
function ^'(r) = e'"^"^. After integrating over the 
angular coordinate (i.e. performing a Hankel transform) 
one obtains 

/•oc 

|$(fc)| / dr r\^{r)\J^{kr) , (34) 
Jo 

where we have taken the modulus to remove constant 
phase factors, and Jv{kr) is the Bessel function of order 
v. As we have already seen, for small v one can approx- 
imate the equilibrium density in position space with the 
TF profile I0J. In this case the integral can be rewritten 
in terms of the variable ( = (r^ - R\ / 2) / {R^_ / 2) as 

(35) 

The limit of this integral for a narrow annulus (r/ ^ 1) 
is of particular interest. A series expansion of the Bessel 
function in powers of rj then gives the simple result 

\m\ = \[^[Mm+o{v^)] , (36) 

where R = l/\/2A is the radius of the annulus. Note 
that odd powers of rj disappear due to parity and that 
the accuracy of this expansion decreases for large k. 




2 4 6 8 10 12 

kR 

FIG. 7: Scaled momentum distributions 128A|$(fc)| V(37r7?) 
as a function of kR (where R — l/\/2A) for (a) v — Q, 
and (b) 1/ = 1 (A = 0.00322, -q = 0.3). Solid lines plot the 
(squared) integral l|84^ calculated with solutions of the Gross- 
Pitaevskii equation, while the dashed lines show Eq. (I35I I 
for the Thomas-Fermi approximation. For comparison, the 
squared Bessel functions [J^(/c7?)]^ are plotted with dotted 
lines. 



Fig. 13 shows the cylindrically-symmetric momentum 
distributions |<I>(fc)p for ^ = and v — 1 and previ- 
ously used parameters A = 0.00322, rj = 0.3. The solid 
line corresponds to a numerical evaluation of Eq. 134|l . 
where the wavefunction is found by solving the time- 
independent Gross-Pitaevskii equation, while the dashed 
line plots the TF result (|35|l . These are scaled for com- 
parison to the Bessel function J^{kR). One sees that 
the Bessel function closely approximates the TF result 
at small momenta, with differences appearing at larger 
kR. Further calculations show that the agreement de- 
grades for larger 77, as would be expected, although it is 
still reasonably good for kR < 5. 

There is also a significant difference between the TF 
and full GP solutions, mostly in the amplitudes of the 
peaks. This is despite the fact that ^ ~ 0.1 d so the con- 
dition for the validity of the TF approximation is well 
satisfied. Crucially, however, the positions of the min- 
ima and maxima of the momentum distribution are rel- 
atively insensitive to the details of the wavefunction. In 
particular, for both ly — and 1^=1 the first few zeros 
and maxima of |<I>(A:)p from the GP solution agree well 
with the roots of the corresponding Bessel functions and 
their derivatives. By virtue of this close correspondence, 
the positions of these peaks and troughs could provide a 
means to measure vorticity in experiments. For example, 
the first maximum could be particularly useful, since its 
position shifts to higher k as v increases. 



VI. SUMMARY 

We have studied a dilute Bose-Einstein condensate 
in a "Mexican hat" trapping potential of quartic minus 
quadratic form. We use the Thomas-Fermi (TF) approx- 
imation to analytically find the equilibrium properties of 
both non-rotating and rotating condensates, and com- 
pare to numerical solutions of the Gross-Pitaevskii (GP) 
equation. In particular, we find that a non-rotating con- 
densate has an annular structure, with a hole in the cen- 
ter, for r] < 1, where 77 = {12gX^/Try^^. We also find 
that for condensates in a frame rotating at sufficiently 
low angular velocities a stable multiply-quantized vortex 
with circulation v > I is present at the center. At higher 
angular velocities our Gross-Pitaevskii simulations show 
that the ground state can instead consist of a centered 
multiply-quantized vortex surrounded by rings of singly- 
quantized vortices, some of which are visible within the 
annulus. 

We have also investigated the dynamics of a conden- 
sate with and without a multiply-quantized vortex, par- 
ticularly the collective modes. A hydrodynamic equa- 
tion was derived within the TF approximation, with so- 
lutions that agree well with the results of evolving the 
GP equation for large interaction strengths and num- 
bers of atoms where the TF approximation is expected 
to be valid. We find that for a non-rotating conden- 
sate in the limit of a thin annulus (ry 0) the mode 
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frequencies are uj = \/ j [j + 1)1 where the number of ra- 
dial nodes j = 0,1,2,... defines famihes of modes con- 
taining different angular quantum numbers m. In this 
limit modes with the same j but different m are degen- 
erate, but diverge as one increases 77. The j — family 
for m > are low-lying modes corresponding to sound 
waves directed around the annulus. Thus a simple an- 
alytical expression for their frequencies can be derived 
from (jj = cq, where c = r// \/l2A is the speed of sound, 
and q — \m.\/R with R = l/\/2A the radius of the annu- 
lus, giving ui = \m\r]/V6. 

The mode frequencies in the presence of a multiply- 
quantized vortex were also analyzed, where it was found 
that the frequency difference between opposite m modes 
in the rotating frame is much smaller than the splitting 
introduced when transforming to the laboratory frame, 
especially for the low-lying modes. The possibility of 
using this splitting for experimental detection of vortic- 
ity is explored, in particular the excitation of low-lying 
m — ±2 modes in a thin annulus. Although in theory 
the splitting could be observable as a "precession" of the 
mode around the annulus, the practicality of this tech- 
nique may be limited due to the difficulty of observing 
the associated small density fluctuations. Hence, as an al- 
ternative method we also investigated the free expansion 
of the condensate after release from the trap. We found 
that the hole in the annulus tends to be filled during the 
expansion, meaning that a vortex is clearly discernible as 
a hole in the expanded density profile. As a consequence 
of the repulsive interactions, the density of the expanded 
cloud is relatively uniform near to the center, allowing 



the width of the vortex core to be used as a measure of 
the vorticity. 

In contrast, when interactions are switched off during 
the expansion we find that the central density is less uni- 
form, and at long times it tends towards the momentum 
distribution of the condensate. Since the momentum dis- 
tribution can be also probed by Bragg spectroscopy, we 
have explored the possibility of using such a measure- 
ment to detect vortex states. We find that for thin an- 
nuli in the Thomas-Fermi approximation the distribution 
takes the form (x [J,y{kR)]'^. Fourier transform- 

ing the stationary solutions of the Gross-Pitaevskii equa- 
tion changes the amplitudes, but not the positions, of the 
peaks, potentially allowing the minima of an experimen- 
tally measured momentum distribution to be compared 
to the zeros of the Bessel function in order to distinguish 
different i' states. 

Finally, we would like to stress that most of the results 
obtained in the thin annulus limit are expected to be 
valid for generic ring-shaped condensates, independent of 
the specific form of the trapping potential. In particular, 
near to the minimum a generic Mexican hat potential 
can be approximated by a harmonic well in the radial 
coordinate, allowing a precise mapping onto the results 
for our quartic minus quadratic potential. 

Acknowledgments 

This work was supported by the Ministero 
dellTstruzione, deU'Universita e della Ricerca. 



[1] M. H. Anderson, J. R. Ensher, M. R. Matthews, 
C. E. Wieman, and E. A. Cornell, Science 269, 198 
(1995). 

[2] C. C. Bradley, C. A. Sackett, J. J. ToUett, and 

R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995). 
[3] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van 

Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, 

Phys. Rev. Lett. 75, 3969 (1995). 
[4] S. Gupta, K. W. Murch, K. L. Moore, T. P. Purdy 

and D. M. Stamper-Kurn, Phys. Rev. Lett. 95, 143201 

(2005). 

[5] V. Bretin, S. Stock, Y. Seurin, and J. Dalibard, 

Phys. Rev. Lett. 92, 050403 (2004). 
[6] S. Stock, V. Bretin, F. Chevy, and J. Dahbard, Euro- 

phys. Lett. 65, 594 (2004). 
[7] A. L. Fetter, Phys. Rev. A 64, 063608 (2001). 
[8] E. Lundh, Phys. Rev. A 65, 043604 (2002). 
[9] K. Kasamatsu, M. Tsubota, and M. Ueda, Phys. Rev. A 

66, 053606 (2002). 
[10] G. M. Kavoulakis and G. Baym, New J. Phys. 5, 51.1 

(2003). 

[11] A. D. Jackson, G. M. Kavoulakis, and E. Lundh, 

Phys. Rev. A 69, 053619 (2004). 
[12] A. L. Fetter, B. Jackson, and S. Stringari, Phys. Rev. A 

71, 013605 (2005). 



[13] J. -K. Kim and A. L. Fetter, Phys. Rev. A 72, 023619 
(2005). 

[14] H. Fu and E. Zaremba, cond-mat/0508515 

[15] M. Cozzini, A.L. Fetter, B. Jackson, and S. Stringari, 

Phys. Rev. Lett. 94, 100402 (2005). 
[16] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 59, 

2990 (1999). 

[17] J. Tempere, J. T. Devreese, and E. R. I. Abraham, 

Phys. Rev. A 64, 023603 (2001). 
[18] E. Nugent, D. McPeake, and J. F. McCann, Phys. Rev. A 

68, 063606 (2003). 
[19] A. Aftalion and I. Danaila, Phys. Rev. A 69, 033608 

(2004). 

[20] P. O. Fedichev and G. V. Shlyapnikov, Phys. Rev. A 60, 
R1779 (1999). 

[21] P. Rosenbusch, V. Bretin, and J. Dalibard, 

Phys. Rev. Lett. 89, 200403 (2002). 
[22] Yu. Kagan, N. V. Prokof'ev, and B. V. Svistunov, 

Phys. Rev. A 61, 045601 (2000). 
[23] A. B. Bhattacherjee, E. Courtade, and E. Arimondo, 

J. Phys. B: At. Mol. Opt. Phys. 37, 4397 (2004). 
[24] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Wahaven, 

Phys. Rev. Lett. 85, 3745 (2000). 
[25] A. Gorlitz, J. M. Vogels, A. E. Leanhardt, C. Raman, 

T. L. Gustavson, J. R. Abo-Shaeer, A. P. Chikkatur, 



11 



S. Gupta, S. Inouye, T. Rosenband, and W. Ketterle, 
Phys. Rev. Lett. 87, 130402 (2001). 

[26] L. P. Pitaevskii and S. Stringari, Bose-Einstein Conden- 
sation (Clarendon Press, Oxford, 2003). 

[27] R. J. Donnelly, Quantized Vortices in Helium II (Cam- 
bridge University Press, Cambridge, 1991). 

[28] F. Zambelli and' S. Stringari, Phys. Rev. Lett. 81, 1754 
(1998). 

[29] F. Chevy, K. W. Madison, and J. Dahbard, 
Phys. Rev. Lett. 85, 2223 (2000). 

[30] P. C. Haljan, L Coddington, P. Engels, and E. A. Cornell, 
Phys. Rev. Lett. 87, 210403 (2001). 

[31] The j = m = solution is excluded since it does not 
conserve particle number, and is therefore unphysical. 

[32] The exact solution of Eq. 1251 can be formally written as 
a series expansion in powers of t], i.e. ujj^ = ]^j, Wjmkif , 
SSjmiC) = e™** J]j,pj>„fc(C)??*, where the functions 
(C) are polynomials of order j + k and parity 
(—1)-'"'"'°. Keeping terms up to ry^, for the most interest- 
ing cases j = 0,1 one has = ni^if /Q and uJ^rn ~ 2 + 
(m'^/2 — 3)?7^/5, with the corresponding phase variations 
SSomiO = [Po(C) + P2(C)m2'7V36]e™'^ and 5Si™(C) = 
[A(C) - P2(C)'7/2 + (mV4 + 6)P3(C)'?V25]e™^ 

[33] W. H. Press, B. P. Flannery, S. A. Teukolsky, and 
W. T. Vetterling, Numerical Recipes in Fortran (Cam- 
bridge University Press, Cambridge, 1992). 



[34] P. Engels, L Coddington, P. C. Haljan, and E. A. CorneU, 
Phys. Rev. Lett. 89, 100403 (2002). 

[35] At the leading order in rj and for small angular veloci- 
ties, where 5njm{C,) ~ Pj(0^^"^'^ ^ the quadrupole opera- 
tor is x^-y^ = [Po(C)+'?A(C)]cos(20)/(2A) = (5no,+2 + 
Sno.~2)/{4:X)+O{ri), mainly exciting the (j — 0,m = ±2) 
modes. The (j = l,m — ±2) modes, whose excita- 
tion amplitude is proportional to r/, are responsible for 
the small-amplitude high-frequency oscillations visible in 
Fig.Hb). 

[36] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dal- 
ibard, Phys. Rev. Lett. 84, 000806 (2000). 

[37] Experimentally, one could realize a purely 2D expansion 
by retaining the confinement in the axial direction. A 
full 3D expansion, inducing a lowering of the effective 
2D coupling constant, would enhance the increase of the 
healing length, changing the scaling of the vortex core 
radius. However, the qualitative behavior, and in partic- 
ular the factor of two difference between the vortex radii 
for V = Q and = 1, is expected to be preserved. 

[38] J. Stenger, S. Inouye, A. P. Chikkatur, D. M. Stamper- 
Kurn, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 
82, 4569 (1999). 

[39] F. Zambelli, L. Pitaevskii, D. M. Stamper-Kurn, and 
S. Stringari, Phys. Rev. A 61, 063608 (2000). 



